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TRAFFIC ANALYSIS IN RANDOM DELAUNAY TESSELLATIONS AND 

OTHER GRAPHS 
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Abstract. In this work we study the degree distribution, the maximum vertex and edge flow in 
non-uniform random Delaunay triangulations when geodesic routing is used. We also investigate 
the vertex and edge flow in Erdos-Renyi random graphs, geometric random graphs, expanders and 
random fe-regular graphs. Moreover we show that adding a random matching to the original graph 
can considerably reduced the maximum vertex flow. 
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Random graphs constitute an important and active research area with numerous applications to 
geometry, percolation theory, information theory, queuing systems and communication networks, to 
mention a few. They also provide analytical means to settle prototypical questions and conjectures 
j> that may be harder to resolve in specific circumstances. In this work we analyze the traffic conges- 

tion on random Delaunay triangulations with non- uniform density, Erdos-Renyi random graphs, 
random fc-regular graphs and random geometric graphs. 

Given a set of points in Euclidean space, the Delaunay triangulation gives a canonical triangulation 
whose vertices are these points. These triangulations have many nice combinatorial and geometric 
properties that make them extremely useful. Moreover, they can also be constructed in Riemannian 
manifolds. However, they do not exist for arbitrary sets of points: certain density requirements are 
needed to ensure that the triangulation can accurately represent both the topology and geometry 
of the manifold. These triangulations are canonically determined by the set of points and have 
many of the properties that Delaunay triangulations have in the Euclidean space. 



The definition of Delaunay tessellations is the same in Riemannian geometry as it is in W l . They 
are defined as having the empty circumscribing sphere property: the minimal radius circumscribing 
sphere for any simplex contains no vertices of the tessellation in its interior. However, there are 
several possible problems with this definition. How do we know the circumscribing sphere is unique? 
For instance, a necessary requirement is that all simplices of our Delaunay tessellation and their 
neighbors are contained inside strongly convex balls. Otherwise, we would run into problems with 
the tessellation being well defined. Since in this work we are not is interested in studying Delaunay 
triangulations per se but in understanding their traffic congestion and graph properties, we focus on 
the following construction. Let M C M? be a either the open unit disk or the whole Euclidean plane 
with density p(x, y), and let V be a Poisson point process of density p on M. The distribution and 
density of these points are determined by the function p. Now let T be the Delaunay triangulation 
of V with respect to the Euclidean metric in R 2 . In this work we are mostly interested in the case 
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where the density is rotationally symmetric (i.e. p(x, y) = f(x 2 + y 2 )) and where 

/ p(x, y)dxdy = oo. 
Jm 

The last condition guarantees that the points in the set V are infinite with probability one. Our 

interest is in understanding how the graph's characteristics change as we change the density. For 

instance, let T be as before and choose xq an arbitrary vertex. Consider the sets 

T n := {x G T : d(x , x) < n} 

with the graph metric (hop metric). For each pair of nodes, consider a unit flow that travels through 
the minimum path between nodes so that the total flow in T n is equal to \T n \ {\T n \ — l)/2- If there is 
more than one minimum path for some pair of nodes v, w, the flow splits among them in a locally 
equal manner: If v' is a vertex on a such a path, the flow out of v' is split equally among neighbors 
{x : d(v' , x) = 1 and d(x, w) = d(v' , w) — 1}, where d is the hop metric. 

Given a node v we define T n {v) as the total flow generated in T n passing through v. In other words, 
T n (v) is the sum off all the geodesic paths in T n which are carrying flow and contain the node v. 
Let My be the maximum vertex flow 



M v (n) := max<T n (v) : v € T n [• 



It is easy to see that for any graph \T n \ — 1 < M n < \T n \(\T n \ — l)/2. Analogously, we can 
define M e (n) as the maximum edge congestion. Note that the same question can be formulated by 
replacing the infinite graph T with a sequence of graphs {G n } c ^' =1 such that \G n \ — > oo. 

One of our main motivations is to understand how M v is affected by changes in the function p. 
We also study how the density affects the degree distribution of the corresponding triangulation. 
Another topic is the maximum vertex flow in Erdos-Renyi random graphs, random geometric graphs 
and random fe-regular graphs. 

It was observed in many complex networks, man-made or natural, that the typical distance between 
the nodes is surprisingly small. More formally, as a function of the number of nodes N, the average 
distance between a node pair typically scales at or below 0(log(N)). Moreover, many of these 
complex networks, specially communication networks, have high congestion. More precisely, there 
exists a small number of nodes called the core where most of the traffic pass through. In this work, 
we present a possible solution to this problem. Given a graph G = (V,E), consider the graph H 
constructed by adding a random maximal matching to G. More precisely, we choose a pair of nodes 
at random from G and add an edge between these two nodes if they are not already connected. 
Now, we remove these two nodes from the possible candidates and repeat this process until there 
is only one or no nodes remaining (depending on the parity of |V|). It is clear that the new graph 
H has the same nodes as before, and it has at most \V\/2 extra edges. Moreover, we added at 
most an extra edge at every node. We show that if the original graph satisfies certain hypotheses 
(exponential growth) then we can reduce the maximum vertex flow to 0(A r poly(log(A r ))). 

This paper is organized as follows. In Section[2j we present the construction of the random Delaunay 
triangulations and we show a necessary condition to guarantee maximum vertex flow of the order 
0(iV 2 ). We also show that for any planar graph, congestion cannot be smaller that ©(A'" 3 ' 2 ) 
regardless of how the flow is routed. In Section [3j we analyze the maximum flow and degree 
distribution for several examples. In Section [4l we analyze maximum vertex and edge flow in 



TRAFFIC ANALYSIS IN RANDOM DELAUNAY TESSELLATIONS AND OTHER GRAPHS 3 

random Erdos-Renyi random graphs, random geometric graphs and random A;-regular graphs. We 
further present the main result on the congestion after a random matching. Finally, Section [5] 
discusses some of the algorithms used, their complexity and how they were implemented. 



2. CONSTRUCTION 

Consider M = B(0,t) C R 2 the open ball of radius t = 1 or t = oo. Let p be a rotationally 
symmetric density on M. Let V be a Poisson process with intensity p on M, and let T be 
the Delaunay triangulation of the set of points in V . In particular, introducing geodesic polar 
coordinates r, 9 makes it apparent that for any Borel set ACM the expected number of points in 
A is 

/ p(r)rdrd6. (2.1) 

J A 

Therefore, if 

2vr / p(r)rdr = oo (2.2) 

Jo 

then almost surely T induces and infinite graph (triangulation) on M. 

Example 2.1. If t = 1 and 

4A 



p{r) 



(1-r 2 ) 2 ' 

then p corresponds to the constant density on the unit ball with the metric of the Hyperbolic space 
under the Poincare disk representation with curvature —A 2 . 

Example 2.2. If t = oo and 

p{r) = 1 
then p corresponds to the constant density on the Euclidean space. 

For notational simplicity, we assume that G V\ otherwise we choose the closest point to the origin. 
Let {r n }'^' =1 be an increasing sequence of positive numbers such that r n < t and bim^oo r n = t. 
Consider T n the graph induced by the finite set {x G T : d(0, x) < r n } where d(x, y) for x,y G T 
is the hop metric in the triangulation. In this Section, we are interested in analyzing the traffic 
behavior on the graph T n as n increases as we discussed in the introduction. For instance, if Gm 
is Kn, the complete graph with N vertices, then M V (N) = N — 1. It was proved, in [5j [6] that 
if G is an infinite Gromov hyperbolic graph then M v {n) = @(iV 2 ). (For more details in Gromov 
hyperbolic graphs, see p2]). In particular, trees have congestion of the order iV 2 . 

One of our main questions is, what is the rate of growth of M v (n) and M e {n) as n — > oo. In 
particular, under what conditions is 

lim *«g > 07 (2.3) 

ra-»oo TV (n) z 

The next Theorem gives a sufficient condition, but first let us recall a result from Riemannian 
geometry. Given a rotationally symmetric Riemannian surface with metric ds 2 = f{r) 2 dr 2 + 
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g(r) d6 , its curvature K{r) is equal to 

K(r)= - 1 i(m) (24) 

Theorem 2.3. Assume t is equal to 1, and let p be a differentiable, strictly increasing function 



that is strongly convex (i.e., hip is convex) and satisfies condition (2.2). Then there exists a > 
such that 

limsup-^ = a>0. (2-5) 

Proof. Let M = B(0, 1) be the Riemannian surface with metric ds 2 = f{r) 2 dr 2 + r 2 f{r) 2 d9 2 . Then 
its area form dA is equal to dA = rf{r) 2 drd9. Hence, the Euclidean un it ba ll -6(0, 1) with volume 



density p(r) = f(r) corresponds to this metric. By applying equation (2.4) to M we see that 



-1 / f'(r) d 2 \ 

Since / is increasing and strongly convex, we observe that K(r) < for all < r < 1. Let V 
be a Poisson process in B(0, 1) with intensity p. This is equivalent to a Poisson process in M of 
constant unit intensity. Let T be its Delaunay triangulation, and let {r n }^ =1 be an increasing 
sequence converging to 1. Denote by T n the restriction of T to B(0,r n ), and let N(n) = \T n \- 

As before, we assume that G V (otherwise we choose the closest point to the origin). For each 
< S < 1, the volume vol(-B(0,5)) < oo and vo\(B(0, 5) c ) = oo. Therefore, we have only a finite 
number of nodes in B(0, 5) and an infinite number outside this ball. Since we are interested in 
the asymptotic behavior of M v (n), we can ignore all the traffic flow generated inside B(0, 5) for n 
sufficiently large. Given x and z in B(0,5) c , let ^ xz be the geodesic in M connecting these two 
points. Since M has negative curvature, the geodesies bend over to the origin. In particular, let 
uj be the number such that vol(5(0, w)) = 1 and let Cy be the shadow of this ball seen through 
the point y to the unit circle according to the Euclidean metric and let C y be the shadow of this 
point according to the metric on M as shown in Figure [TJ Since M has negative curvature, it is 
clear that Cy C C y . Moreover, the set C y depends only on the distance from y to the origin. Let 
v(r) = length(Cj / ) in the unit circle dM = S . 

By construction, the points in V are distributed in such a way that there is one point per unit area 
in M. Therefore, it is clear that the geodesies in T are quasi-geodesics in M. Then the traffic flow 
in T n is equal to the traffic in (B(0, r n ), ds) up to a constant. Hence, up to a fixed constant, 

limsuplimsup — -. > limsuplimsup " > 2ir / v(r)dr > (2-6) 

0<<5<1 n->oo N(ny <<5<1 n^oo N(n) z J s 

proving our claim. □ 

As mentioned above, the rate of growth of M v (n) is between ®{N) and Q(N 2 ) for any graph with 
A^ nodes. The next Theorem shows that this function can be narrowed down if the graph is planar. 

Theorem 2.4. If G is a planar graph, then 

e(N 3/2 ) < M v < e(N 2 ) (2.7) 
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Figure 1. Geodesic path and shadows in M. The shadow Cfj is in magenta and 
the shadow C„ is in red. 



where N is the number of nodes in G. Moreover, the same result holds for any traffic routing (not 
necessarily geodesic routing). 



Proof. By the planar separator theorem (PQ) in any N- vertex planar graph G = (V, E), there exists 
a partition of the vertices of G into three sets A, S and B such that each of A and B has at most 
2N/3 vertices, S has 0{yN) vertices, and there are no edges with one endpoint in A and one 
endpoint in B. It is not required that A or B form connected subgraphs of G. The set S is called 
the separator for this partition. Therefore, all the traffic between nodes in A and B has to pass 
through a node in S. Hence, by the pigeonhole principle there exists a node in S with traffic at 
least 9(iV 3 / 2 ). □ 



It was proved in j2] that if the graph G has exponential growth then 

®{N 2 /\og(N)) < M v < B{N 2 ) 
independently of how the traffic is routed. 



2.1. Poisson Process Construction. In this subsection, we describe a way to construct a real 
ization of the Poisson V in practice. Let 

1)2 . TTT)2 



$ : 



2.9) 



be the measure preserving map given in polar coordinates by <&(re ld ) = a(r)e ie where a : [0, oo) — > 
[0, t) is an increasing function that depends on the function p. This map is measure preserving if 
and only if for every x > 0, 

ra(x) 



vol R 2(B(0,x)) = TTX^ 



2tt / p{y)ydy =: 2irF(a(x)). 
/o 
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Figure 2. Delaunay triangulation with A equal to 1 (left), 5 (center) and 10 (right). 



Therefore, a has to be equal to 

a(x) = F^~ 1 \x 2 /2) (2.10) 

where F^ 1 ' is the inverse with respect to composition of the function F. Let 7-r2 be a Poisson point 
process on the Euclidean plane of constant intensity 1. We can construct V as the set of points 
^T-j^). Hence, everything boils down to constructing a Poisson process on the Euclidean space 
which can be done by partitioning the space into a set of disjoint annuli and uniformly distributing 
points in these annuli. 



3. Congestion Analysis and Degree Distribution 

In this Section, we show some simulation results obtained for the maximum edge and vertex flow 
as well as the degree distribution for different densities p. 



3.1. Hy perb olic Plane with the Poincare Disk Representation. Assume that p is as in 

2A and F(x) 



Example 
hence 



2.1 



Then it is an easy calculation to show that J p{y)y dy 

( „2 \ 1/2 

a(x) 



i-ir" 



2Xx z 
1-x 2 ' 



4A + x 2 I ' 

In Figure [2j we show the Delaunay triangulation with 10,000 nodes of the previous surface for 
different values of the intensity A. In table [T] we show the diameter, average and maximum vertex 
and edge flow. We observe that as A increases, the maximum vertex and edge flows decrease but 
the averages increase. 

Table 1. Traffic Congestion 



A 


Diameter 


AvgVflow 


MaxVflow 


AvgEflow 


MaxEflow 


1 


20 


130764 


3.57E+07 


22031 


8.64E+06 


5 


31 


190195 


2.94E+07 


31906.8 


6.56E+06 


10 


37 


221776 


2.50E+07 


37145.9 


7.53E+06 
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Figure 3. Delaunay triangulation for N = 10,000 and a = 0.01 (left). Case 
Hyperbolic space representation in M 2 (center) and Euclidean case (right). 

3.2. Hyperbolic Space Representation in R 2 . We can put a density p in the entire R 2 that 
corresponds to the hyperbolic metric on R 2 . More specifically, consider R 2 with the metric 

ds 2 = dr 2 + sinh (^) M 2 

that corresponds to the hyperbolic space with uniform curvature —K. Hence, if K = 1, we get the 
classical hyperbolic space. The area form is equal to dA = sinh(r) drd6. Then taking p{r) = ^ r ' 
we get a density that corresponds to the density to the Hyperbolic space density. In this case, 

F(x) = cosh(x) — 1 

and it is easy to see that 

a{x) = acosh(x 2 /2 + 1). 

3.3. A "Generalized" Hyperbolic Space. Let t = 1 and p(r) = 2r t ° ° where a > 0. Note 
that a = 1/2 corresponds to the Poincare disk. It can be show that f p(y)ydy = 2 1/(1 and 

2x i/a 



F(x) 
hence 

a(x) = 



1 - xV a 



x 2 



A + x 2 



Figure [3] gives a realization for N = 10, 000 and a = 0.01 in comparison with the Euclidean density 
and the density corresponding to the Hyperbolic space representation in M 2 . Figure |4J shows that the 
degree distribution of its corresponding Delaunay triangulation follows a power law for a = 0.01 
and 0.05. Moreover, the degree of the power law distribution can be tuned up by choosing the 
appropriate intensity. Therefore, this is a way to construct a random triangulation with a power 
law degree distribution. Not surprisingly, a = 0.5, which correspond to the Poincare disk, does not 
follows a power law distribution. 
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Degree Distribution "Generalized" Hyperbolic 




Figure 4. Degree distribution for N = 10, 000 and a = 0.01, a = 0.05 and a = 0.5. 




Figure 5. Delaunay triangulation with intensity 0.1 (left), intensity 1 (center) and 
intensity 10 (right). 

3.4. Another Example. Let t = 1 and p(r) = ,._ y S for A > 0. Then it can be show that 



A 2 (r-l/2) 



JT=f^ 



f p(r)rdy = ( \_J 2 and F(x) = \ (1 Z x)2 , hence 



(1-X)2 



a(x) 



\ + x' 



Figure [5] shows the Delaunay triangulation of this for different values of the intensity A. 
3.5. Slow Growing. Let t = oo and p(r) = r +1 \ - Then 

F(x)= p(y)ydy = ln(x + l). 
Jo 



Hence, 



a(x) = exp I — - 1 I . 
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Maximum Geodesic Vertex Flow 



Maximum Geodesic Edge Flow 
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Figure 6. Log-log plots of the maximum vertex congestion for different densities 
as a function of the number of nodes, averaged over 20 trials. 



In Figure [6j we see the asymptotic growth of the maximum vertex flow and the maximum edge flow 
as a function of the number of nodes N of T n for the densities from Subsections 3.1-3.3 It can be 



shown that by changing the function p, the maximum vertex conges tion changes from any possible 
growth rate between @(iV 3 ' 2 ) and @(N 2 ). As discussed in Theorem 2.4 these are the only allowed 
growths for planar graphs. 



4. Erdos-Renyi Graphs, Random Geometric Graphs and Expanders 



In this Section, we explore the maximum vertex congestion with geodesic routing for Erdos-Renyi, 
random geometric graphs and expanders. More specifically, let ER(iV,pjv) be the random Erdos- 
Renyi graph with N nodes and p^ = 2log(N)/N. Note that these graphs are connected almost 
surely since the probability pn is slightly above the connectivity threshold (1 + e) \og{N)/N . Let 
RGG(A r , rjy) be the random geometric graph in the unit square with n nodes and radius r^ = 
y i /2\og(N)/N. These graphs again are connected since the radius is slightly above the connectivity 
radius. Now we turn to a particular construction of random expander graphs. We will assume that 
k = 21 is an even integer. To build a fc-regular graph on iV vertices {1, . . . , iV}, what we do is pick I 
permutations m, . . . , ix\ : {1, . . . , N} — >■ {1, . . . , N}, and let Qk,N be the graph formed by connecting 
v to 7Tj(i>) for all v e {1, ... , N} and i = 1, . . . , 1. Formally, this is not always a A;-regular graph 
since there could be multiple edges and loops. It is well known that eliminating these problematic 
edges the graph is a fc-regular and an expander with probability tending to 1 as N — > oo. 

Figure [7[ shows the maximum vertex flow and the maximum edge flow as a function of the number 
of nodes for these three graph families. In the case of the expander family, we use k = 6. As can 
be appreciated from these plots, there is very little congestion even for geodesic routing. Moreover, 
note that at the connectivity threshold is the worst case scenario since if we enlarge p^ or rjy the 
congestion decreases even more. We should also point out that the expander graphs have the least 
edge and vertex congestion even though these graphs have the least number of edges among the 
three families. 

As discussed in Section [2l Gromov hyperbolic graphs have congestion of the order 0(iV 2 ). In 
particular, any fc-regular tree (also called a Bethe lattice) has highly congested nodes. More 
precisely, the following result was proved in [5]. 
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Maximum Geodesic Vertex Flow 
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Figure 7. Log- log plot of the maximum vertex and edge congestion as a function 
of the number of nodes and averaged over 20 realizations. 

Proposition 4.1. For a k-regular tree with N nodes, 



M v 



k-\ 
2k 



(N - l) 2 + N - 1. 



(4.1) 



Furthermore, the maximum congestion occurs at the root. 



Note that trees are some of the most congested graphs one can consider. The reason is that much 
of the traffic must pass through the root of the tree. This leads to the following question. 

Question 4.2. Is there a way to reduce the vertex and edge congestion in a graph by adding a 
relatively small number of edges? 



We answer this question positively, but we first need the following Lemma. 

Lemma 4.3. Let G be a graph with bounded geometry, i.e. sup„ eG deg(t> ) < A < oo. Then for 
every v £ G, the traffic flow passing through v satisfies 

T{v) <A 2 {A-l) D - 2 -D 2 

where D = diam(G). 



(4.2) 



Proof. Let v £ G and define Sk := {x £ G : d(v, x) = k). Then it is clear that G = {v} U lL=i Sp 
and moreover |Sfc| < A(A — l) fe ~ 1 . Furthermore, 

T(v)< £ |S fc |.|5,| 

k+l<D 

where the inequality is coming from the fact that if k + I > D then the geodesic path between a 
node in S& and a node in 5; does not pass through v. Hence, 



T(v)< ^2 A 2 {A-l) k+l - 2 < A 2 (A-L 



D ~ 2 D 2 . 



k+KD 



□ 



Given a graph G, a matching M in G is a set of pairwise non-adjacent edges; that is, no two edges 
share a common vertex. The following result is due to Bollobas and Chung |8j . 
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Theorem 4.4. Suppose G is a graph on N vertices with bounded degree k satisfying the property 
that for any x £ V{G), the i-th neighborhood Ni(x) of x (i.e., Ni(x) = {y ■ d(x,y) = i}) contains 
at least ck(k — l) l ~ 2 vertices for i < (1/2 + e) log fc _ 1 (iV) J where e and c are fixed positive values. 
Then by adding a random matching to G gives a graph H whose diameter D(H) satisfies 

log k (N) -C< D(H) < \og k {N) + log, log(iV) + C 

with probability tending to 1 as N goes to infinity, where C is a constant depending on e and c. 

Now we are ready to answer the previous question. 

Theorem 4.5. Let G be a graph on N vertices with bounded degree k satisfying the hypothesis 
of Theorem \4-4\ Let H be the graph constructed by adding a random matching to G. Then the 
maximum vertex flow M V {N) with geodesic routing on H satisfies 

M V {N) < k c Nlog(N)logl(N) + o(Nlog(N)log 2 k (N)). (4.3) 

with probability tending to 1 as N —> oo. 



Proof. Applying Lemma 4.3 gives 

M V (H) < k D ^D{H) 2 . 

Since 

£>(#)< log fe (iV) + log fe log(Af) + C 

with high probability, we see that 

M V (H) < (log fe (iV) + log fc log(iV) + C)2 A .log fc (JV)+log fc log(7V)+C 

= k c Nlog{N)log 2 k (N) + o{Nlog(N)log 2 k (N)). 

D 

In particular, in every /c-regular tree, we can reduce the congestion from @(N 2 ) to 0(iVpolylog(iV)) 
by adding at most an extra edge for every node. The same applies for the Delaunay triangulations 
constructed in the previous sections. 

It was proved in |7J that the diameter of a random /c-regular satisfies the following property. 

Theorem 4.6. For k > 3 and sufficiently large N, a random k-regular graph G with N nodes has 
diameter at most 

log fc _ 1 (iV) + log fc _ 1 (log fe _ 1 (iV)) + C 
where C is a fixed constant depending on k and independent on N. 

Hence by reasoning as in the proof of Theorem |4.6| we obtain the following corollary. 



Theorem 4.7. Let G be a random k-regular graph. Then the maximum vertex flow M V {N) with 
geodesic routing on G is smaller than 

k c N\og\{N) + o{Nlogl(N)). (4.4) 

with probability tending to 1 as N —> oo. 
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Maximum Geodesic Vertex Flow 



Maximum Geodesic Edge Flow 
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— BL 
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f(N) = Nlog(N) 




Figure 8. Edge and vertex flow for a Bethe lattice of degree 6 as a function of the 
number of nodes, and the same after adding one or two random matchings. 

In Figure pi we see the maximum vertex and edge flow for a random regular tree of degree 6 (Bethe 
Lattice of degree 6) as a function of the number of nodes, and the same after adding one or two 
random matchings to the edge list of this graph. By adding one random matching to the edge list, 
we reduce the flow from order N 2 to Nlog(N) 2 . 

5. Implementation and Algorithms 



Section [2?l] explains that we use a Poisson process to generate random points and remap the distance 
from the origin via r — > a{r). Then we need the Delaunay triangulation with a flow computation 
based on routing over geodesic paths. Thus the main computation can be broken down into four 
steps: 

(1) Generating appropriately distributed pseudorandom vertices, 

(2) Computing the Delaunay triangulation, 

(3) Finding the shortest path between each pair of points, 

(4) Computing the traffic flow. 

The purpose of this section is to describe these steps, consider their run time and space complexity, 
and study the efficient implementation of the critical the ones that consume the most time and 
space. 

Steps [T] and [2] can be used alone to get degree statistics for graphs that are too big for the traffic 
flow computation, so it is noteworthy that Step [I] handles an n-vertex graph in time 0(n), while 
Step [2] requires 0(n log n). 

Step [I] simply selects pseudorandom points (xi, yi) uniformly in a disk, converts to polar coordinates, 
replaces r by a(r), and converts back to rectangular coordinates. If the Poisson rate suggests n 
points, we pseudorandomly select a number N from the rate n Poisson distribution, and select N 
points in the disk. For maximum flexibility, our C++ implementation allows the user to enter the 
function a(r) in postfix notation using arithmetic operations and standard transcendental functions. 
This adds a constant factor overhead to the run time for Step [TJ but this is subsumed by Step [2] 

The Delaunay triangulation for Step [2] is computed via Fortune's sweep line algorithm running as 
a separate process [IT] . 
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Step [3] could potentially use a lot of time and space, so efficiency is important. The all-points 
shortest path algorithm just takes the graph's adjacency matrix A and computes A p for a power p 
not less than the graph diameter using min-plus arithmetic For our application, A is very sparse 
and p is only O(logn), so there is no harm in exponentiating via repeated multiplication by A. 
(The process stops when such a multiplication has no effect.) 

The matrix for the all-points shortest path problem keeps track of shortest known path lengths, 
and sparseness means having a lot of infinite entries. As suggested by Figure [9j it is important for 
the min-plus matrix multiplication to take advantage of sparsity, yet be efficient when the matrix 
is dense. Hence, we store it as a dense (integer) matrix, but also have a bitmap that keeps track 
of sparsity. A big switch statement can consider 8 entries at a time, and do 0-8 add and minimize 
operations as needed. Furthermore, one 32-bit boolean operation can cause 32 entries to be skipped 
if the matrix is very sparse. It also helps to compute the load factor and switch to dense-matrix 
arithmetic when it gets high. 




Figure 9. The matrix load factor at each step in the all-points shortest path com- 
putation for an n = 800 problem with diameter 14. 

Let D be the matrix that gives the path lengths from the all-points shortest path computation. 
The shortest-path routing that forms the basis for Step [4] requires flow from a vertex v destined for 
a vertex w to go to neighbors Vi of v for which D nw = D vw — 1. When there is more than one such 
neighbor, we divide the flow into equal parts. Hence, each edge from v has a sparse vector that 
tells, for each possible destination, what fraction of v's outbound traffic takes that edge. These 
redistribution vectors are sparse, and their non-zero entries are of the form \/k for small positive 
integers k. Table [2] tells how these k values are typically distributed. 

Table 2. Popularity of k values in redistribution vectors 



1 


2 


3 


4 


5 


6 


7 


> 8 


49% 


33% 


12% 


3.9% 


1.0% 


0.57% 


0.13% 


0% 



Table [2] also provides an estimate of the sparsity in a typical redistribution vector. Since 1 • 0.49 + 
2 • 0.33 + . . . 7 • 0.0013 = 1.76 and the tabulated data are for a planar graph with average degree 6, 
the load factor is 1.76/6 = 0.293. Hence one byte almost always suffices to give a non-zero entry 
and specify how many zeros since the last non-zero. 
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The traffic flow computation for Step HI begins with a starting flow matrix S such that each Si 7 j 
gives the initial specification for flow from vertex % to vertex j. The object is to compute an overall 
flow matrix 

F = S + ujS + uj 2 S + uj 3 S + --- , (5.1) 

where a; is a linear operator formed by componentwise multiplication of columns of a flow matrix 
times redistribution vectors, e.g., the redistribution vector r^- for an edge i,j adds diag(rjj) times 
column i of w's argument matrix to column j of w's result matrix. Note that ( |5.1[ ) specifies a fixed 
point of the iteration F <— S + ujF. This fixed point is reached after O(logn) steps because io v is 
the zero operator whenever exponent p exceeds the graph diameter. 

To compute some column j of S 4- cjF, we need redistribution vectors r^ for every i £ Ej, where 
Ej is the set of all i such that i,j is in the edge set E: 



S J + 5Z & a &(. r iiJ) F i- 



This is trivial, except that we can reduce the need for scratch memory by copying the result back 
to Fj as soon as the old Fj is no longer needed. 

Now consider the total memory requirement for the flow computation on a graph of n vertices and 
m edges. The all-points shortest path matrix D (including sparsity bits) requires (4 + g)n 2 bytes. 
Since i,j and j, i require different redistribution vectors at about 0.3n bytes each, the total memory 
for redistribution vectors is about 0.6mn bytes. The starting flow S requires little or no memory. 
Finally, F and the associated scratch memory require between 8n s and 16n 2 bytes. For a planar 
graph with m ~ 3n, the total should between 13. 8n 2 and 21. 8re 2 , and experiments have shown that 
19n 2 is typical for such planar graphs. 

The run time for the all-points shortest path computation depends on sparsity. The adjacency 
matrix A has 1m non-zeros, and the load factors in Figure [9] sum to 7.5, a number that does not 
seem likely to grow much as n increases. Hence, the all-points shortest path computation requires 
about 15mn add-and-minimize operations. Each F <— S + ojF computation requires one multiply- 
add for each non-zero in the redistribution vectors, for a total of about 0.6mn. Finally, the 0.6mn 
multiply-add's are repeated d + 1 times, where the graph diameter d is O(logn). 
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